
clear

SaveCon = 'Results/full_CF';
load([SaveCon '.mat']);
SaveRes = 'Results/full.mat';
load(SaveRes);

% Recover FID_Grid
optCount = opt; optCount.SelSubSample = 0; 
optCount.TrunkAge = 0; % Option to trucate age (added later, so may not be present in some results files)
[FID_Grid,Action,Age,AgeState,States,pproc,DS] = prepdataM10(optCount);

% Make sure FID_Grid match:
GridNMatch = any(CF.FID_Grid - FID_Grid ~= 0);
if GridNMatch
    disp(' CF grid does not match with estimates option generated grid ')
end


%% Elasticities decomposition (Article Table 3: Long-run elasticities)

% Individual acreage and yield elasticities:
IndAcreageElas = ((CF.LandUseCF05./CF.LandUseBase) - 1)/0.05;
IndYieldElas   = (((CF.OutputDetailCF05./CF.LandUseCF05)./(CF.OutputDetailBase./CF.LandUseBase))-1)/0.05;
IndOutputElas  = ((CF.OutputDetailCF05./CF.OutputDetailBase) - 1)/0.05;

TableElas.AcreageSweight = nansum(IndAcreageElas.*CF.OutputDetailBase)/nansum(CF.OutputDetailBase);
TableElas.AcreageAweight = nansum(IndAcreageElas.*CF.LandUseBase)/nansum(CF.LandUseBase);
TableElas.YieldSweight   = nansum(IndYieldElas.*CF.OutputDetailBase)/nansum(CF.OutputDetailBase);
TableElas.Output         = nansum(IndOutputElas.*CF.OutputDetailBase)/nansum(CF.OutputDetailBase);

TableElas = struct2table(TableElas)

%% Graph short-medium run elasticities

Time = 1:150;
for t = Time
    Output(t)  = (mean(CF.OutputPathCF05(t,:),2)/mean(CF.OutputPathBase(t,:),2) - 1)/0.05;
    Yield(t)   = (mean(CF.YieldPathCF05(t,:),2)/mean(CF.YieldPathBase(t,:),2) - 1)/0.05;
    Acreage(t) = (mean(CF.AcreagePathCF05(t,:),2)/mean(CF.AcreagePathBase(t,:),2) - 1)/0.05;
end

figure1 = figure(1); % Short and long run
% Create axes
axes1 = axes('Parent',figure1,'XTick',[1 3 5 7 9 11 13 15 20 25 30 40 50 60 70 80 90]);
xlim(axes1,[1 90]);
box(axes1,'on');
grid(axes1,'on');
hold(axes1,'all');

% Create multiple lines using matrix input to plot
plot1 = plot([Output(1:90)' Yield(1:90)' Acreage(1:90)'],'Parent',axes1,'LineWidth',3);
set(plot1(1),'DisplayName','Output');
set(plot1(2),'DisplayName','Yield');
set(plot1(3),'DisplayName','Acreage');

% Create xlabel
xlabel('years');

% Create ylabel
ylabel('elasticity');

% Create legend
legend1 = legend(axes1,'show');
saveas(figure1, [SaveCon '_longrun'], 'epsc')

figure2 = figure(2); % Short run
% Create axes
axes1 = axes('Parent',figure2,'XTick',[1 3 5 7 9 11 13 15 17]);
xlim(axes1,[1 17]);
box(axes1,'on');
grid(axes1,'on');
hold(axes1,'all');

% Create multiple lines using matrix input to plot
plot2 = plot([Output(1:17)' Yield(1:17)' Acreage(1:17)'],'Parent',axes1,'LineWidth',3);
set(plot2(1),'DisplayName','Output','Color','black','LineStyle','-');
set(plot2(2),'DisplayName','Yield','Color','black','LineStyle','--');
set(plot2(3),'DisplayName','Acreage','Color','black','LineStyle',':');

% Create xlabel
xlabel('years');

% Create ylabel
ylabel('elasticity');

% Create legend
legend2 = legend(axes1,'show');
saveas(figure2, [SaveCon '_shortrun'], 'epsc')

% Land use classification:
LU = nan(size(DS,1),1);
LU(DS.mb_uso_detalhe == 3 | DS.mb_uso_detalhe == 9)  = 1; % Forests
LU(DS.mb_uso_detalhe == 4 | DS.mb_uso_detalhe == 12) = 2; % Savanna (including grasslands)
LU(DS.mb_uso_detalhe == 15)  = 3; % Pasture
LU(DS.mb_uso_detalhe >= 19 & DS.mb_uso_detalhe <= 21) = 4; % Agriculture
LUclasses = unique(LU(~isnan(LU)));
k = length(LUclasses);
TableExpansionDecomposition.LU             = LUclasses;
TableExpansionDecomposition.Area           = zeros(k,1);
TableExpansionDecomposition.ShareArea      = zeros(k,1);
TableExpansionDecomposition.Elasticity     = zeros(k,1);
TableExpansionDecomposition.TotElasContrib = zeros(k,1);
TableExpansionDecomposition.ShareOfExpan   = zeros(k,1);


u = 1;
for w = LUclasses'
    TableExpansionDecomposition.Area(u)           = sum(LU == u);
    TableExpansionDecomposition.ShareArea(u)      = sum(LU == u)/sum(~isnan(LU));
    TableExpansionDecomposition.Elasticity(u)     = nansum(IndAcreageElas(LU == u).*CF.LandUseBase(LU == u))/nansum(CF.LandUseBase(LU == u));
    TableExpansionDecomposition.TotElasContrib(u) = nansum(IndAcreageElas(LU == u).*CF.LandUseBase(LU == u))/nansum(CF.LandUseBase);
    TableExpansionDecomposition.ShareOfExpan(u)   = TableExpansionDecomposition.TotElasContrib(u)/TableElas.AcreageAweight;    
    u = u + 1;
end

TableExpansionDecomposition = struct2table(TableExpansionDecomposition);
TableExpansionDecomposition.Properties.RowNames = {'Forests', 'Savanna', 'Pasture', 'Cropland'};
TableExpansionDecomposition




